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Collagen is the main structural and load-bearing 
element of various connective tissues, where it 
forms the extracellular matrix that supports cells. 
It has long been known that collagenous tissues 
exhibit a highly nonlinear stress-strain relation¬ 
ship Q], Hj], although the origins of this nonlinear¬ 
ity remain unknown [3]. Here, we show that the 
nonlinear stiffening of reconstituted type I colla¬ 
gen networks is controlled by the applied stress, 
and that the network stiffness becomes surpris¬ 
ingly insensitive to network concentration. We 
demonstrate how a simple model for networks 
of elastic fibers can quantitatively account for 
the mechanics of reconstituted collagen networks. 
Our model points to the important role of nor¬ 
mal stresses in determining the nonlinear shear 
elastic response, which can explain the approxi¬ 
mate exponential relationship between stress and 
strain reported for collagenous tissues [l|. This 
further suggests new principles for the design of 
synthetic fiber networks with collagen-like prop¬ 
erties, as well as a mechanism for the control of 
the mechanics of such networks. 

Significance 

We report nonlinear rheology experiments on collagen type 
I networks, which demonstrate a surprising concentration in¬ 
dependence of the network stiffness in the nonlinear elastic 
regime. We develop a model that can account for this, as well 
as the classical observations of an approximate exponential 
stress-strain relationship in collagenous tissues, for which a 
microscopic model has been lacking. Our model also demon¬ 
strates the importance of normal stresses in controlling the 
nonlinear mechanics of fiber networks. 
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Introduction 

Collagen type I is the most abundant protein in mam¬ 
mals where it serves as the primary component of many 
load-bearing tissues, including skin, ligaments, tendons 
and bone. Networks of collagen-type I fibers exhibit 
mechanical properties that are unmatched by man-made 
materials. A hallmark of collagen and collagenous tissues 
is a dramatic increase in stiffness when strained. Qualita¬ 
tively, this property of strain stiffening is shared by many 
other biopolymers, including intracellular cytoskeletal 
networks of actin and intermediate filaments HHZ|. On 
closer inspection, however, collagen stands out from the 
rest: it has been shown that collagenous tissues exhibit a 
regime in which the stress is approximately exponential 
in the applied strain |1]. The origins of this nonlinearity 
are still not known [3|,|gj, and existing models for biopoly¬ 
mer networks cannot account quantitatively for collagen. 
In particular, it is unknown if the nonlinear mechanical 
response of collagen originates at the level of the individ¬ 
ual fibers HUES® or arises from nonaffine network 
deformations as suggested by numerical simulations HU¬ 
ES- 

Here, we present both experimental results on reconsti¬ 
tuted collagen networks, as well as a model that quanti¬ 
tatively captures the observed nonlinear mechanics. Our 
model is a minimal one, of random networks of elastic 
fibers possessing only bending and stretching elasticity. 
This model can account for our striking experimental 
observation that the stiffness of collagen becomes inde¬ 
pendent of protein concentration in the nonlinear elas¬ 
tic regime, over a range of concentrations and applied 
shear stress. Our model highlights the importance of 
local network geometry in determining the strain thresh¬ 
old for the onset of nonlinear mechanics, which can ac¬ 
count for the concentration independence of this thresh¬ 
old that is observed for collagen 1E3, in strong con¬ 
trast to other biopolymer networks. Finally, our model 
points to the important role of normal stresses in deter¬ 
mining the nonlinear shear elastic response, including the 
approximate exponential relationship between stress and 
strain reported for collagenous tissues [1]. 

Results and Discussion 

In contrast to most synthetic polymer materials, 
biopolymer gels are known to exhibit a strong stiffening 
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response to applied shear stress, in some cases leading to 
a more than 100 -fold increase in the shear modulus, at 
strains as low as 10 % or less, before network failure 0 - 
7]. Here, we perform rheology experiments on reconsti¬ 
tuted networks of collagen type I, a key component of 
many tissues. We measure the differential shear modu¬ 
lus K = dcr/d^f relating the shear stress a to the strain 7 . 
We plot this in Fig. [l^i as a function of the applied stress. 
At low stress (and strain), we observe a linear elastic re¬ 
sponse with K = G, the linear shear modulus, which 
increases with collagen concentration. These networks 
also exhibit a strong increase in their stiffness K above 
a threshold stress that increases with concentration. Re¬ 
markably, for network concentrations ranging from 0.45 
to 3.6 mg/ml, the modulus becomes insensitive to con¬ 
centration in the nonlinear regime, where K increases 
approximately linearly with a: here, for a given sam¬ 
ple preparation (e.g., polymerization temperature), the 
various K vs a curves overlap, in spite of the fact that 
the linear moduli of these samples vary by two orders of 
magnitude. 

Moreover, the approximate linear dependence of K on 
a in our reconstituted networks is consistent with the 
empirically established exponential dependence of stress 
on strain in collagenous tissues [l|, since cr oc exp ( 7 / 70 ) 
implies that K = da /dy oc a. While qualitatively similar 
stiffening with applied stress has been reported for other 
biopolymers MSGIGl, both the linear dependence 
of K on a and the insensitivity of the nonlinear stiffening 
to network concentration appear to be unique to collagen. 


Physical picture. Although surprising at first sight, the 
features seen in Fig.[]Ji can be understood in simple phys¬ 
ical terms for athermal networks of fibers that are soft 
to bending and where the nonlinear network response is 
controlled by stress. At low stress, if the elastic energy is 
dominated by soft bending modes, the linear shear modu¬ 
lus G should be proportional to the fiber bending rigidity 
ft. Of course, G also depends on the density of collagen, 
as can be seen in Fig. [lji. The concentration can be char¬ 
acterized in geometric terms by p, the total length of fiber 
per unit volume. Since ft has units of energy x length, 
while G has units of energy per volume, we expect that 
G oc ftp 2 [ 2 (| El. Since stress has the same units as 
G, similar arguments apply to the characteristic stress 
ao oc ftp 2 , above which the response becomes nonlin¬ 
ear. For ft = 0, such networks become entirely floppy 
and their rigidity depends on other stabilizing effects, or 
fields , including applied stress [19, I 22 I l25j . Thus, when 
the applied stress a becomes large enough to dominate 
the initial stability due to fiber bending resistance, it is 
expected that K will increase proportional to cr, in a way 
analogous to the linear dependence of magnetization on 
field in a paramagnetic phase. Combining these obser¬ 
vations, one obtains an approximate stiffening given by 
K oc G x (<r/(To) m , where the slope m = 1 . Here, since 
G and ao have the same dependence on concentration, 
one obtains a nonlinear stiffness K that becomes insensi¬ 


tive to concentration. Interestingly, this behavior is nei¬ 
ther expected nor observed for F-actin and intermediate 
filament networks, which are not bend-dominated and 
exhibit a stronger nonlinear stiffening regime, in which 
K oc cr 3 / 2 0,0. 


Model. To test this simple physical picture, as well as 
uncover the mechanisms of collagen elasticity in more 
detail, we study simple/minimal computational models 
of fiber networks, specifically, 2D and 3D lattice-based 
networks [26l428| and 2D Mikado networks pj], H, ' 29j |. 
It is known that the mechanical stability and rigidity 
of networks depends on their connectivity, which can 
be characterized by the coordination number z, defined 
by the number of fiber segments meeting at a junction. 
Prior imaging of collagen networks [30] report an aver¬ 
age connectivity z ~ 3.4. Importantly, this places such 
networks well below the isostatic or critical connectivity 
of 2 = 4 in 2D or z = 6 in 3D required for mechanical 
stability of networks with only spring-like stretching en¬ 
ergies [3l|. As a result, the linear elastic properties are 
expected to be governed by other energies, such as fiber 
bending [ll|, [l2|,[l6|, llOl 12 llj . a s well as by the distance of 2 
from its critical value [22|,[25|. Thus, we generate our net¬ 
works within a range of 2 , straddling the experimentally 
relevant values. Specifically, our 2D and 3D lattice-based 
networks are created with z = 3.2 and our Mikado net¬ 
works have z = 3.6. As we show below, properties such 
as the linear modulus G and the strain threshold for the 
onset of nonlinear elasticity depend on z, although the 
overall form of the nonlinear regime is unaffected. 

In our model, as in our experiments, we impose a 
volume-preserving simple shear strain 7 and minimize 
the total elastic energy H of the network, consisting of 
the sum of elastic energies of the individual fibers. The 
elastic energy of a fiber is calculated using a discrete form 
of the extensible wormlike-chain model that accounts for 
both local stretching and bending [29] (also Supporting 
Information). The network stiffness K is calculated as 


K = 


1 d 2 n 

V <9 7 2 5 


(i) 


where V is the volume of the system. Since K depends 
on the energy per unit volume, and the energy involves 
an integral along the contour of all fibers in the system, 
K is naturally proportional to the total length of fiber 
per volume, p, which is proportional to the protein con¬ 
centration c. Thus, K can be expressed as (Supporting 
Information) 


K = upK. ( 7 , k) , (2) 

where p is the fiber stretching modulus and k = k,//j,£q 
is a dimensionless measure of the relative bend-stretch 
stiffness, with £q a measure of the spacing between fil¬ 
aments. Here, p oc £l~ d . For lattice-based networks, 
we define £q to be the lattice spacing, while for Mikado 
networks we use the average distance between crosslinks. 
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Fig. 1. Stiffening of reconstituted collagen type t networks, (a) Differential shear modulus K vs shear stress a for reconstituted 
collagen type t networks at varying protein concentrations and different polymerization temperatures (red: 37 °C, blue: 25 °C). 
Lines of unit slope serve as visual guides. The filled red diamonds represent a network polymerized at 37° C at a concentration 
of 1.8 mg/ml with 0.2% glutaraldehyde cross-linkers. The inset is a schematic of a typical stress vs strain curve indicating the 
stiffness K as the tangent or differential shear modulus and the point ( 70 , 00 ) at the onset of stiffening, (b) Simulation results 
for 2D (red) and 3D (green) lattice-based and 2D Mikado (blue) networks for various reduced bending rigidities k = 10 -2 (■), 
k — 4 x 10 -3 (o), k — 2 x 10 -3 (•), k = 10 -3 (v), k = 6 x 10 -4 (♦), k — 2 x 10 -4 (0), and k — 10 -4 (▼). The lattice-based 
networks (red and green) have connectivity z = 3.2 corresponding to an aspect ratio L/£q = 5, while the Mikado networks 
(blue) have z = 3.6 with L/£o = 11. We see that changing z affects the overall magnitude of the moduli, but not the functional 
form of the stiffening response. The inset shows stiffness vs stress curves from a 3D lattice-based network simulation under 
volume-preserving extension, where T is the extensional stress and A is the extension ratio, (c) Comparison of K vs a curves 
obtained from experiment (a) and 3D lattice-based network simulation (x) under shear. Multiplicative factors for the stiffness 
and stress axes have been chosen for coincidence of the linear modulus and the stress at the onset of nonlinearity, (d) A 3D 
confocal image of a reconstituted collagen type I network shows a highly branched local geometry (right). Collagen fibers are 
hierarchically assembled of fibrils (diameter: lOnm) which in turn consist of staggered collagen molecules (diameter: 1.5nm). 
The overall fiber diameter is of order lOOnm, which makes the fibers sufficiently rigid enough to be modeled as an elastic beam, 
(e) Confocal images show differences in network geometry at different polymerization temperatures. Polymerizing collagen at 
25 °C creates networks of straighter, less branched fibers in contrast to networks polymerized at 37 °C. (f) The 2 D network 
geometries used in the simulations. 
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The shear stress a can be expressed in a similar fashion 
as a = fipTt (7, ft). For a given network structure, JC and 
E are dimensionless functions of only 7 and k. 

In our simulations, we determine both K and a for 
various values of ft. We do this for networks with p = 1 
and £q = 1 . Thus, our simulation values of both moduli 
and stress are in units of p/i ^ 1 in d dimensions. We plot 
K vs a in Fig. [lb>. For an elastic rod of diameter 2 a and 
Young’s modulus E, the parameter k is proportional to 
the fiber volume fraction (ft, since ft = na A E/A, p = tt a 2 E 
and k = a 2 /( 4 £q) oc (ft [ll], Eli. We thus consider 
values around k < 10 -3 to compare with experiments, 
where the protein volume fraction varies over a range of 
approximately one decade around 0 . 1 %. 

Consistent with our experiments, our model networks 
also show an approximately linear relationship between 
stiffness K and shear stress a, as shown in Fig. [Tb> [ 26 | . 
We also study networks under extension, for which our 
model predicts a linear relationship between the stiffness 
and extensional stress, as shown in the inset to Fig. []Jd. 
Thus, our model can also account for prior experiments 
on collagenous tissues, which report such a linear rela¬ 
tionship [ 1 ]. Moreover, both experiments and theory 
show a very surprising result in the stiffening regime , 
where the if vs a curves for different networks are seen 
to cluster around a common line, and where networks of 
varying protein concentrations exhibit the same stiffness 
at a given level of applied shear stress; i.e., the network 
stiffness K becomes independent of network concentra¬ 
tion and appears to be governed only by the applied stress 
in the nonlinear regime. 

For low stress, the linear regime is indicated by a con¬ 
stant stiffness K = G, for which our model predicts the 
linear dependence on k: G oc pk oc ftp 2 . This is consis¬ 
tent with both our observed increase of G with collagen 
concentration in the experiments (Fig. 153b). as well as 
with prior reports showing an approximate quadratic de¬ 
pendence of G on concentration BE- Moreover, to 
test whether for a given concentration G increases with 
ft, we show data with glutaraldehyde (GA) cross-linkers, 
which increases the bending rigidity of collagen fibers [ 32 ] 
(Fig. UK). Not only are these results consistent with the 
predicted increase in G, but the K vs a curve still collapse 
onto the corresponding data for non-GA cross-linked net¬ 
works in the stiffening regime. Thus, our model can ac¬ 
count for the features observed in the experiments. For 
a more direct comparison we plot theoretical and exper¬ 
imental stiffening curves together in Fig. [TJ3. Moreover, 
both 2 D and 3 D results exhibit similar behavior, suggest¬ 
ing that stiffening is independent of dimensionality for a 
given local network geometry (Fig. [Tb>). 

In the nonlinear regime, the observed independence of 
K/a on concentration, and therefore on the typical spac¬ 
ing between fibers, suggests that the stiffening should 
be understood purely in geometrical terms, and quanti¬ 
ties such as the characteristic strain 70 at the onset of 
stiffening should be independent of sample parameters 
such as concentration and ft. Figure [ 2 b shows that 70 
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Fig. 2. Independence of characteristic strain 7 at the onset 
of stiffening on concentration. (Color online) (a) Onset strain 
70 obtained from simulations vs fiber bending rigidity k. (b) 
Experiments showing independence of 70 on protein concen¬ 
tration. (c) The upper panel shows the stiffness vs strain in 
a 2D lattice over the broad range of ft, and reveals a strain- 
stiffening regime highlighted by the shaded region. The inset 
shows the same data from Fig. [lb normalized by concentra¬ 
tion and plotted vs strain, with the dashed lines correspond¬ 
ing to the 70 values and the shaded regions highlighting the 
stiffening regimes. In both simulation and experiment, 70 
is estimated as the strain at which the stiffness is roughly 
twice the linear modulus. The lower panel shows the overall 
contribution of bending energy to the total elastic energy in 
the network. (d,e) Experimental data for 37 °C showing the 
strain dependence of the raw stress and stiffness. The symbols 
denote the same concentrations as shown in Fig. Ob. 


is indeed independent of ft, and is thus independent of 
both p and ft, throughout the range k < 10 -3 . The 
strain threshold 70 does, however, depend on the the 
connectivity, z, as well as the type of the network, i.e., 
whether lattice-based or Mikado. As we show in the Sup¬ 
porting Information, in the strongly bending-dominated 
limit, our model predicts a simple scaling dependence of 
70 oc ( £o/L ) 2 on the aspect ratio L/£ 0, where L is the 
average length of the fibers. In general, 70 decreases 
with increasing L/t 0 or z. For a given network type, 
lattice-based or Mikado, this aspect ratio is an entirely 
equivalent measure of connectivity to z: there is a one- 
to-one relationship between these two quantities, which 
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increase (decrease) together. By construction, our net¬ 
works have local coordination numbers strictly less than 
4, which also represents the physiological upper bound 
of two fibers crossing at a cross-link. As the aspect ratio 
L/£ 0 oc, £ ^ 4 corresponding to the limit of very long 
fibers cross-linked many times to each other. Conversely, 
as 2 decreases toward 3 (a branched structrure) the as¬ 
pect ratio decreases toward unity. Thus, stiffening in our 
model networks is controlled by geometry, specifically via 
the aspect ratio L /£o or equivalently, the coordination 
number z. 

Collagen is known to form branched network struc- 
tures (3(1 |33j (see Fig. [TJi), whose pore size scales as 
1 / y/c [34[ . Changing the concentration only changes the 
degree of branching while preserving the local geometry, 
including the aspect ratio; i.e., networks at different con¬ 
centrations look alike, apart from an overall scale factor. 
The onset strain 70 is then predicted to be independent 
of concentration, and indeed we observe this experimen¬ 
tally (Fig. [2b). Although this is consistent with prior 
experiments on collagen ||, [Tt} , it is in strong contrast to 
reports for other biopolymer networks [4 IM Tl9l |35| . 

Role of local network geometry. We can now un¬ 
derstand quantitatively the features in our experiments 
based on three key assumptions: (i) the networks are 
athermal, (ii) are bend-dominated and (iii) their geom¬ 
etry at different concentrations is self-similar, i.e., the 
network structures at different concentrations are scale- 
invariant in that they are characterized by the same (as¬ 
pect) ratio L/£o- We test the last hypothesis by prepar¬ 
ing collagen networks with different geometries. The 
structure of collagen networks strongly depends on the 
polymerization conditions, such as pH, ionic strength or 
temperature [ 9 I, IhgI HJot ] (see Fig. [lb). Changing the local 
geometry, and specifically L/t 0 by changing the poly¬ 
merization temperature does not affect the form of the 
stiffening response nor the collapse of the data in the 
nonlinear elastic regime, in either the model or the ex¬ 
periments, apart from a change in the ratio K/cr. The 
stiffening curves of networks with different geometries 
cluster around distinctly different curves of approximate 
unit slope (Fig. [lb). Moreover, less branched networks 
show a lower 70 (Fig. [2b)- This is consistent with simu¬ 
lation results when comparing Mikado with lattice-based 
networks (Fig. [TJd and Fig. [2b). To confirm that this is 
due to network geometry, and not to the temperature 
at which the rheology measurements are performed, we 
polymerize a network at 37 °C and subsequently cool it 
to 25 °C. We then perform rheology measurements at 
25 °C and find that despite its larger linear modulus, the 
stiffening regime coincides with networks polymerized at 
37 °C, demonstrating that network geometry, indeed, sets 
the prefactor K/cr (Fig. IS3bh 

To understand the stiffening mechanism, we first exam¬ 
ine which of the two modes, stretching or bending, dom¬ 
inates the stiffening regime. Prior work has suggested 
that stiffening corresponds to a transition from bending- 


to stretching-dominated behavior 0 . Our simulations 
show that bending is dominant throughout the stiffening 
regime (Fig. [2b, Fig. |S5]) . When stretching modes finally 
become dominant, all /C vs 7 curves converge, as shown in 
Fig. [2b. In most cases, this only occurs after the network 
stiffness has increased by more than an order of magni¬ 
tude. Moreover, when stretching dominates, we find a 
distinct stiffening behavior characterized by if ^ cr 1 / 2 
(Fig. |S 6 |) [Hj. Thus, we find three distinct rheological 
regimes: ( 1 ) a linear elastic regime, ( 2 ) a bend-dominated 
stiffening regime, and (3) a stretch-dominated stiffening 
regime. Interestingly, the approximate K ~ a regime we 
observe in our collagen networks is consistent with the 
second of these regimes, which occurs before the transi¬ 
tion from bend- to stretch-dominated behavior. 

The existence of a distinct bend-dominated nonlinear 
regime and the corresponding concentration-independent 
nonlinear response in Figs, [lb and b depends crucially 
on the sub-isostatic nature of the networks, as well as 
on small values < 10“ 2 of k in the model and volume 
fraction <p in experiments. The collagen networks we 
study here are, indeed, all sub-isostatic with respect to 
stretching alone [HI, [3l|, since z ~ 3.4 lies well below 
the critical connectivity of 6 in 3D (4 in 2D) at which 
pure spring networks first become stable. As either k 
or the aspect ratio L/t 0 increase, a transition to stretch- 
dominated linear elastic behavior is expected, even in 3D, 
where the networks remain clearly sub-isostatic 00 . 
However, over the range 2.5 < L/t 0 < 11 that we study 
here (Figs, [lb and S4), which includes reported collagen 
network structures, we consistently see effects of bend- 
dominated network response in our model, including the 
concentration-independent nonlinear behavior. Here, k 
acts as a stabilizing interaction or field for networks in 
their linear elastic regime, with G oc ft. The intermedi¬ 
ate nonlinear regime, where we find K ~ a in our simu¬ 
lations and experiments, can be understood in terms of 
marginal stability together with the stabilizing effect of 
applied stress. 


Normal stresses. Biopolymer networks, including col¬ 
lagen, have been shown to develop large negative normal 
stresses [2f| 0, |4l| . This is in contrast to most elastic 
materials that exhibit positive normal stresses, as first 
demonstrated by Poynting [42[, who showed elongation 
of wires under torsion. Biopolymer gels have been shown 
to do the opposite. In experiments, the constraint nor¬ 
mal to the sample boundaries leads to the build-up of 
tensile stress at these boundaries when simple shear is im¬ 
posed. These normal stresses can stabilize sub-marginal 
networks. In Fig. [3b, we show that the linear modulus 
grows in direct proportion to an applied normal stress. 
We hypothesize that the network stiffness could arise 
from the normal stresses that develop under shear strain 
due to the imposed constraint at the boundaries: 

K ~ Gq + X a N 


(3) 
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Fig. 3. Stiffening induced by normal stresses. (Color online) 
(a) The change in the linear shear modulus grows in direct 
proportion to an external normal stress cfn applied on the 
shear boundaries, (b) Comparison of K (black) with Gq+x&n 
(blue) vs shear stress in the linear and stiffening regimes. The 
local slope a = 1 (red dashed line) in the stiffening region is 
shown as a visual guide. The upper insets show the variation 
of ct as a function of bending rigidity/protein concentration. 
The lower inset shows the reduction in network stiffness when 
the normal stress is removed. 


where Go is the linear shear modulus in the absence of 
any normal stress o/v and y is a constant. In Fig. [3t>, 
we show a direct comparison of K and Go + x a N vs 
cr, where dy is independently measured in our simula¬ 
tions. The linear regime is characterized by Go in the 
absence of a at- In the stiffening regime, there is excellent 
agreement between K and Go + x a Ni and both show the 
same local slope a ~ 1 consistent with the unit slope 
in Figs, [lji and [TJd. Finally, we confirm our hypothesis 
by performing further relaxation of the networks when 
the normal stresses are released. In the lower inset of 
Fig. [3 }d, by removing the normal stresses, we observe a 
significant reduction of the stiffness throughout the stiff¬ 
ening regime. This clearly supports the hypothesis that 
normal stresses control the stiffening of fiber networks 
under simple shear. Moreover, upon closer examination 
of the model predictions, we see a small but systematic 
evolution of the stiffening exponent a with k, as shown 
in the upper left inset of Fig. [SJd. A similar evolution is 
seen in our experiments as a function of concentration, 
as shown in the right panel of the upper inset of Fig. [3t>. 
This agreement between experiment and model further 
justifies our identification of k with network concentra¬ 
tion. 


Concluding Remarks 

The development of normal stresses in these networks 
is intimately related to the volume-preserving nature of 
simple shear deformations, both in our rheology experi¬ 
ments and in our simulations. In our model, removal of 
the normal stress leads to a reduction in the volume of 
the system and a concomitant reduction in the stiffness. 
While collagenous tissues in vivo are subject to more 
complex deformations, approximate volume conservation 
is valid in many cases, e.g., due to embedded cells [l[. 
Our results suggest that any volume-preserving deforma¬ 
tion should lead to similar behavior in stiffness vs stress. 
In particular, in addition to accounting for the approx¬ 
imate exponential stress-strain relationship known em¬ 
pirically for collagen under extension [l|, our model also 
predicts that the nonlinear (Young’s) modulus should be¬ 
come concentration independent for a given extensional 
stress, in a way similar to the case of simple shear. This 
can be seen in the inset to Fig. [lj) and Fig. S7. 

The concentration-independence and collapse of the 
stress-stiffness curves seen in Fig. []Ji appears to be 
unique to collagen networks, at least among biopolymers. 
Within our model, this property depends on three as¬ 
pects: (1) the athermal and simple elastic response of 
the constituent fibers, (2) the bend-dominated response 
of the network in its linear elastic regime, and (3) the 
linear scaling of stiffness with stress, given by K ^ a 171 , 
where m = 1. These properties are also expected for 
collagen type II, another fiber-forming type of collagen. 
For intracellular biopolymer networks, however, the lat¬ 
ter property (3) is strongly violated: for actin and in¬ 
termediate filament networks, a stronger stiffening, with 
m ~ 3/2 is observed. Interestingly, no such collapse or 
concentration-independence has been reported for those 
systems. One interesting consequence of the approxi¬ 
mate collapse of the stress-stiffness relationship, com¬ 
bined with the lack of concentration dependence of strain 
(in Fig. [2ji, b), is that the local deformation of such matri¬ 
ces is expected to become nearly uniform in the nonlinear 
elastic regime, even in the presence of large local inhomo¬ 
geneities in network density. This can have a stabilizing 
effect under excessive mechanical loading. The present 
work has identified the key properties that can form the 
basis for design of biomimetic networks with similar non¬ 
linear properties to collagen. 

The importance of the nonlinear stiffness of collagen 
matrices comes in part from the inherent stability that 
such stiffening can impart to whole tissues: if collagen 
network elasticity were linear, then such networks would 
either fail or have to strain by more than 200-300% un¬ 
der the maximum stresses in our experiments. Moreover, 
an initial soft elastic response of collagen also seems to 
be important physiologically: high stiffness due to exces¬ 
sive collagen production, e.g., during fibrosis, scar forma¬ 
tion or around tumors is known to contribute strongly to 
pathological processes at the cellular level, where it can 
drive the so-called epithelial-to-mesenchymal transition 
and affect cell differentiation. Thus, both a soft initial 
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linear response, as well as a strongly non-linear stiffening 
regime of collagen matrices are important for individual 
cells and tissues. Apart from tissues with high content 
collagen, such as tendon and skin, most soft tissues have 
collagen content in the range of tenths of % to a few %, 
which corresponds to a range from our densest reconsti¬ 
tuted networks up to about a decade higher in concen¬ 
tration [43||. In such tissues, we expect nonlinear effects 
such as we report here to appear at shear stresses of order 
kPa, which is a level of stress easily reached, for instance, 
by traction forces of fibroblasts Q- Thus, we expect the 
kind of stiffening reported here to be relevant to many 
soft tissues. While collagen networks have been known to 
exhibit nonlinear mechanics that is qualitatively similar 
to other biopolymer networks, it has become increasingly 
clear that the underlying mechanism of collagen stiffen¬ 
ing differs from that of other biopolymers Not only 
does the present work shed light on the origins of colla¬ 
gen matrix mechanics, but it can also form a basis for 
the design of synthetic networks to mimic collagen and 
other extracellular matrices for tissue engineering. 

Materials and Methods 

Polymerization of collagen networks. We dilute type 
I collagen (BD Biosciences, San Jose, CA) at 4°C to 
the desired final concentrations between of 0.45 mg/ml 
and 3.6 mg/ml in lx DMEM (Sigma Aldrich, St. Louis, 
MO) with 25 mM HEPES added and adjust the pH to 
9.5 by addition of 1M NaOH. We fill the solution into 
the rheometer geometry preheated to 25 °C or 37 °C as 
indicated and allow for at least two hours of polymer¬ 
ization. To stiffen some samples, we pipette a solution 
of 0.2% glutaraldehyde (Sigma) in lx PBS (Lonza, Al¬ 
lendale, NJ) around the rheometer geometry once the 
networks have polymerized for 45 minutes and incubate 
these samples for three hours before performing experi¬ 
ments. 

Rheometry and data analysis. We perform the exper¬ 
iments on an AR-G2 rheometer (set to strain-controlled 
mode) or an ARES-G2 strain-controlled rheometer (both 
TA instruments, New Castle, DE) both fitted with a 25 
mm PMMA disc as top plate and a 35 mm petri dish 
as bottom plate and set a gap of 400 fi m. We prevent 
evaporation by sealing the samples with mineral oil, ex¬ 
cept for experiments on crosslinked collagen; here, we 
use a custom-built solvent trap, which allows for the ad¬ 
dition of the crosslinking solution. We monitor the poly¬ 
merization of all samples by continuous oscillations with 
a strain amplitude of 0.005 at a frequency of 1 rad/s. 
Subsequently, we impose a strain ramp with a rate of 
0.01/sec and measure the resulting stresses. We fit each 
stress-strain data set with a cubic spline interpolation 
and calculate its local derivative, which we then plot ver¬ 
sus stress. 

Generation of disordered phantom networks. We 

take a W x W triangular lattice or a W x W x W face- 


centered cubic (FCC) lattice of spacing to generate 
the disordered phantom network in two or three dimen¬ 
sions, respectively. In d— dimensions, the lattice occu¬ 
pies a volume V = voW d , where vq is the volume of a 
unit cell. Periodic boundaries are imposed to eliminate 
edge effects. A continuous chain of lattice bonds along 
a straight line forms a single fiber. The lattice vertices, 
having 6-/12-fold connectivity (i.e., coordination num¬ 
ber) in 2D/3D, are freely-hinged cross-links, where fibers 
rotate about each other with no resistance. We reduce 
the average connectivity using the following procedure. 
In a 2D triangular lattice, we randomly select two out of 
the three fiber strands at a vertex on which we form a 
binary cross-link, i.e., with 4-fold connectivity. The re¬ 
maining strand crosses this vertex as a phantom and does 
not interact with the other two strands. This is done at 
every vertex until all cross-links are binary. We further 
dilute the lattice by randomly removing bonds with prob¬ 
ability q = 1 — p, where p is the probability of an existing 
bond. After dilution, fibers that span the system size 
may still be present and these could lead to unphysical 
contributions on the macroscopic network stiffness. To 
avoid this, we make sure that every fiber has at least one 
diluted bond. All remaining dangling ends are further 
removed. Finally, nodes are introduced at the midpoint 
of every lattice bond so that the first bending mode on 
each bond is represented. The procedure just described 
effectively reduces both the average connectivity to z < 4 
and the average fiber length to L = £o/q and generates 
a disordered phantom triangular lattice [26]. A similar 
procedure as described can be implemented on the FCC 
lattice to generate a three-dimensional equivalent [27j . 

Generation of Mikado networks. We generate these 
networks [29] by random deposition of monodisperse 
fibers in the form of rods of length L W onto a 
two-dimensional W x W box, which occupies a volume 
V = W 2 . Each rod’s center of mass (x cm ,p cm ) and ori¬ 
entation ip relative to a fixed axis are each drawn from 
a uniform distribution. The box has periodic boundaries 
such that if a rod intersects any side of the box, it crosses 
over to the opposite side. A cross-link is assigned to the 
point wherever a given pair of rods intersect. Every time 
a rod is deposited, the cross-linking density L/l c is up¬ 
dated, where l c is the average distance between neighbor¬ 
ing cross-links. Deposition stops as soon as the desired 
cross-linking density is achieved, after which all dangling 
ends are removed. Midpoint nodes are introduced on the 
rod between a pair of cross-links. 

Discrete extensible wormlike chain model. The in¬ 
ternal degrees of freedom in the network is the set of 
spatial coordinates {r*} of all discrete nodes (i.e., cross¬ 
links, phantom nodes and midpoint nodes) on every fiber. 
Each fiber in the network is semiflexible, i.e., the elastic 
response to a given deformation is determined by both 
its stretching modulus fi and bending rigidity k. When 
the network is deformed, the nodes undergo a displace- 
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ment {77} —»• {r'}. The extension of a fiber segment 
(ij) between nodes i and j along a fiber is given by 
Mij = t'ij-tij, where 4- = ||r' -r'|| and4,- = ||rj — r-*|| is 
the rest length of the strand. Note that for lattice-based 
networks, iij reduces to the bond rest length £q for all 
(ij). The total stretching energy of a fiber is then cal¬ 
culated by summing up the contributions of a chain of 
strands along its backbone: 


U 


stretch 




<b‘> 



The bending of a fiber segment involves a triplet of con¬ 
secutive nodes (ijk) along the backbone. The local curva¬ 
ture at node j is estimated as ||d£/ds|| ~ Stj = \\ijk~Uj ||, 
where Uj is a unit vector oriented along (ij). The net 
contribution of consecutive segments (ijk) along a fiber 
leads to its bending energy 

"Hbend = J2 Ij (jf) » 

(ijk) \ i J 

where /' = \ (% + Zjk)- Adding up ?4tretch + "Hbend over 
all fibers in the network yields the total elastic energy. 

Rheology simulation. We simulate rheology on the 
networks by imposing an affine simple shear strain 7. We 
fix the fiber stretching modulus fi = 1 and inter-filament 
spacing £ 0 = 1. We vary k to probe a range of bending 
rigidities. We steadily increase the strain in dy steps to 
cover a strain range of 0.1% to 1000%. At each strain 
step, the total elastic energy is minimized by relaxing 
the internal degrees of freedom using a conjugate gradi¬ 
ent minimization technique [ifij. Lees-Edwards bound¬ 
ary conditions are used when calculating the lengths of 
strands that cross the shear boundaries [46j |. From the 
minimum energy 77, we extract the shear stress a and 
differential shear modulus K: 


^ _ 1 on ts — do _ 1 d 2 n 
0 ~ v <97 > A _ y a 7 2 • 
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Supporting Information 

Dimensionless Shear Modulus and Bending 
Rigidity For a homogeneous elastic rod [47] of ra¬ 
dius a and Young’s modulus E , the stretching modulus 
/i = 7T a 2 E and bending rigidity k = ja 4 E. The bending 
length scale defined lul as £\> = yjnjii is a length of order 
the rod diameter, since k/ fi oc a 2 . So for every fiber seg¬ 
ment of length £q, we can express the bending rigidity in 
dimensionless form as k = W = k/jll£ o- The differen¬ 
tial shear modulus is derived from the energy density of 
the network, which requires calculating the total elastic 
energy per unit volume U = y . For a network occupying 
a volume V and composed of N fibers, this can be eval¬ 
uated as a sum of the elastic energies of all semiflexible 
fibers /: 



where d£(s)/ds and \dt(s)/ds\ are the local length change 
and local curvature at a point s on the fiber with unit 
tangent t(s). In a discrete network construction where 
the fibers are divided into a total of n segments (ij), this 
can be approximated as 


N 




£q 


/=i 
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The quantity (£ (7, k))f is a dimensionless elastic energy 
averaged over all fibers in the network, and (£ (7, k)) s is 
averaged over all fiber segments, of which there are n. 
Thus, since n £0 is the total length of fiber in the system, 

U « np(S (7, k)) s , 

where p is the network concentration in total fiber length 
per volume. Differentiating with respect to 7 yields the 
shear stress a = /ipE(7,£), where the dimensionless 
stress 

s ( 7 ,k) = ( 7 ,k)) s . 

Similarly, the dimensionless shear modulus /C = ^ is 
related to the shear modulus by K = pplC (7, k). 

The concentration p is also related to the fiber rigid¬ 
ity k. For any given network structure of stiff rods, 
a segment of length £q and cross-section a 2 occupies a 
volume fraction p oc a 2 £ o/£ 3 > where the typical mesh 
size £ ^ £$. It follows that the concentration of fiber 
material p ~ p oc o? j £\, and since the fiber rigidity 
k = nj p£\ a 2 /4 we obtain k oc p. 



Fig. SI. (Color online) Schematic showing two interacting 
fiber strands fi (blue) and fj (red) as well as other strands 
(gray). Circles denote points of mechanical constraints in the 
form of cross-links or branch points. Relaxation of fi along 
its backbone tugs fj inducing a transverse displacement 8 £± 
(blue arrow) and longitudinal displacements 8£\\ (red arrows). 
In a similar manner, fi experiences both displacements from 
its interaction with other strands. The zoom-in shows a sim¬ 
ple first approximation geometric relation between these dis¬ 
placements. 


the linear elastic response of the network is dominated 
by fiber bending interactions. The backbone relaxation 
of fi induces on fj a transverse displacement 8 £± oc 7 L 
and a longitudinal displacement 8 £\\, as shown on the 
schematic. The longitudinal displacement which is a lo¬ 
cal retraction of /j, is related to the transverse displace¬ 
ment as 


(4-«u) +6f ± =f 0 
S £|[ = 4 - 


Sh 


yi, K 

' s^o ' 


24 

stj 
4 ’ 


such that for small strains, we have S£\\ 


Since 


on average there are L/£ 0 fibers attached to any given 
fiber, the total longitudinal displacement 8 L\\ resulting 
from the backbone relaxations of these other fibers can 
be expressed as 8 L\\ = 8 £\\ ~ ^7^- I n the low strain 

limit 8 L\\ <C 7T, i.e., the total longitudinal displacement 
of the fibers is negligible in comparison to their own back¬ 
bone relaxations. The critical strain 70 is obtained when 
8 L\\ ~ 7T, i.e., 


7 = 7 o ~ l ~j~ 


(Sl) 


which corresponds to the onset of stiffening. Thus, onset 
70 of stiffening is set by the geometric length scale aspect 
ratio L/£ 0 . 

We emphasize that Eqn. [Si] applies to the asymptotic 
bend-dominated limit. This we observe in Fig. [S 2 ] for 



Geometric Dependence of the Critical Strain The 

schematic in Fig. ISTl shows two fiber strands fi and /j, 
each of length £q intersecting at a cross-link. The aver¬ 
age length of the fibers in the network is L. Each fiber 
undergoes a backbone relaxation 7 L and we assume that 


Fig. S2. Critical strain 70 as a function of L/£ 0 . Networks 
with bend-dominated linear elasticity show a shift in the crit¬ 
ical strain with the aspect ratio L/£o described by a line of 
best fit that agrees well with the model. For larger aspect 
ratios, the scaling crosses over to a weaker dependence. 
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L /£o < 5, which is the relevant parameter range in our 
comparison of simulation and experiment. The scaling 
crosses over to a weaker dependence for much larger L/t o, 
where it appears to show L~ l dependence. Such scaling 
has been reported in previous work [27]. However, in con¬ 


trast to the geometric mechanism presented in the cur¬ 
rent work, their L~ l dependence is based on an energetic 
crossover from bending to stretching regimes. Whether 
the scaling changes from L~ 2 to L~ x needs to be further 
investigated. 
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Fig. S3. Concentration dependence of the linear shear modulus and stiffness vs stress curves. (Color online) 
(a) Linear shear modulus vs protein concentration at different polymerization temperatures 37 °C and 25 °C. (b) Linear 

modulus obtained from simulations on networks with different geometries: 2D/3D lattice and 2D Mikado, (c) Stiffness vs stress 
curves normalized by the concentration of collagen networks polymerized at 37 °C. For the network at a concentration of 1.8 
mg/mL, 0.2% glutaraldehyde (GA) cross-linkers are added (filled symbols) to increase the bending rigidity of the fibers, (d) 
Dimensionless stiffening curves from simulations on a 3D lattice for various fiber rigidities, (e) Stiffness vs stress curves showing 
the effect of running the rheology at 25 °C for a network polymerized at 37 °C (solid blue trace). For comparison, we show the 
result when the rheology is run at the same temperature as the polymerization at 37°C (solid red trace). The inset shows the 
increase in linear modulus (black trace) of a network polymerized at 37 °C as the temperature cools down to 25 °C with time 
(blue trace). 
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Fig. S4. Shift of stiffening onset with network geometry. (Color online) Shear stiffening of a 2D phantom network with 
average connectivity z — 3.6 (blue) and z — 3.0 (green). The aspect ratios are L/£o = 5.2 and L/i o = 2.5, respectively, (a). 
Stiffness versus strain shows the shift of the onset of stiffening to a lower value with increasing z or L/i o indicated by the 
arrows, (b) Stiffness versus stress shows an increase in the amplitude ratio K/cr as indicated by the upward shift of the curves 
in the stiffening regime with increasing z or L/i o. The solid line of unit slope serves as a visual guide. 



a/k 


Fig. S5. Stretching and bending contributions to the total energy. The ratio of stretching energy to bending energy 
is less than unity in the stiffening regime, i.e., the shaded region from the critical strain 70 to the strain indicated by the 
thick dashed line. This shows that stretching modes are subdominant to bending in this regime. The insets show the relative 
contributions of stretching and bending to the total elastic energy. 























14 



Fig. S6. Stretch-dominated stiffening. (Color online) Network simulation showing shear stiffening curves for various bending 
rigidities including the zero limit (red dashed curve). This limit corresponds to a network governed purely by stretching modes 
and as the figure shows, it leads to a different stiffening behavior where the modulus scales as . The small deviation from 
the slope of 1/2 at low stress for the k = 0 limit is due to a finite-size effect. The line of unit slope only serves as guide to the 
eye. 



T 


Fig. S7. Stiffening under simple extension. (Color online) Network simulation of stiffening under volume-preserving 
extension, (a) Stiffness vs tensile stress curve for a 3D network with different fiber rigidities. The line of unit slope is shown 
as a guide to the eye. The inset shows a schematic tensile stress vs tensile strain curve and how the stiffness is obtained, (b) 
Stiffening curve on a 2D network under extension with and without volume constraint. 






















